Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod
# Load libraries, install if needed
library(tidyverse); theme_set(theme_light(base_size = 11))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(qwraps2)
library(wesanderson)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB)
# To load entire cache in interactive r session, do:
# qwraps2::lazyload_cache_dir(path = "R/analysis/cpue_model_cache/html")
# Specify map ranges
ymin = 52; ymax = 60; xmin = 10; xmax = 24
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
strip.background = element_rect(fill = "grey95")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'gray10', margin = margin()),
strip.background = element_rect(fill = "gray95"),
legend.position = c(0.7, 0.02),
legend.direction = "horizontal"
)
}
# Make default base map plot
xmin2 <- 303000; xmax2 <- 916000; xrange <- xmax2 - xmin2
ymin2 <- 5980000; ymax2 <- 6450000; yrange <- ymax2 - ymin2
plot_map_raster <-
ggplot(swe_coast_proj) +
xlim(xmin2, xmax2) +
ylim(ymin2, ymax2) +
labs(x = "Longitude", y = "Latitude") +
geom_sf(size = 0.3) +
theme_plot()
plot_map_raster_labels <-
plot_map_raster +
annotate("text", label = "Sweden", x = xmin2 + 0.25*xrange, y = ymin2 + 0.8*yrange, color = "black", size = 1.9) +
annotate("text", label = "Denmark", x = xmin2 + 0.029*xrange, y = ymin2 + 0.43*yrange, color = "black", size = 1.9, angle = 75) +
annotate("text", label = "Germany", x = xmin2 + 0.07*xrange, y = ymin2 + 0.025*yrange, color = "black", size = 1.9) +
annotate("text", label = "Poland", x = xmin2 + 0.55*xrange, y = ymin2 + 0.1*yrange, color = "black", size = 1.9) +
annotate("text", label = "Russia", x = xmin2 + 0.95*xrange, y = ymin2 + 0.2*yrange, color = "black", size = 1.9) +
annotate("text", label = "Lithuania", x = xmin2 + 1*xrange, y = ymin2 + 0.47*yrange, color = "black", size = 1.9, angle = 75) +
annotate("text", label = "Latvia", x = xmin2 + 0.99*xrange, y = ymin2 + 0.68*yrange, color = "black", size = 1.9, angle = 75)
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/mdat_cpue_q_1_4.csv")
# Calculate standardized variables
d <- d %>%
filter(year < 2020 & year > 1992 & quarter == 4) %>%
mutate(oxy_sc = oxy,
temp_sc = temp,
depth_sc = depth,
) %>%
mutate_at(c("oxy_sc", "temp_sc", "depth_sc"),
~(scale(.) %>% as.vector)) %>%
mutate(year = as.integer(year)) %>%
drop_na(oxy, depth, temp)
plot_map_raster +
geom_point(data = d, aes(X*1000, Y*1000))
pred_grid1 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(1_2).csv")
pred_grid <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid_(2_2).csv")
pred_grid <- bind_rows(pred_grid1, pred_grid)
# Standardize data with respect to data grid:
pred_grid <- pred_grid %>%
mutate(year = as.integer(year)) %>%
filter(year %in% c(unique(d$year))) %>%
mutate(depth_sc = (depth - mean(d$depth))/sd(d$depth),
temp_sc = (temp - mean(d$temp))/sd(d$temp),
oxy_sc = (oxy - mean(d$oxy))/sd(d$oxy)) %>% # Need to scale these to the mean and sd in the data!
drop_na(oxy, depth, temp)
spde <- make_mesh(d, xy_cols = c("X", "Y"),
n_knots = 200,
type = "kmeans", seed = 42)
# Plot and save spde
png(file = "figures/supp/density/spde.png", units = "in", width = 6.5, height = 6.5, res = 300)
plot(spde)
dev.off()
# Depth spline + oxy spline
# Takes about 30 minutes
m <- sdmTMB(density ~ 0 + as.factor(year) + s(depth_sc) + s(oxy_sc) + s(temp_sc),
data = d,
mesh = spde,
family = tweedie(link = "log"),
spatiotemporal = "AR1",
spatial = "on",
time = "year",
reml = FALSE,
control = sdmTMBcontrol(newton_steps = 1))
#> Warning: The model may not have converged. Maximum final gradient:
#> 0.0163339061734842.
tidy(m, conf.int = TRUE)
sanity(m)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> x `thetaf` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> x `ln_phi` gradient > 0.001
#> ℹ See `?run_extra_optimization()`
#> ℹ Or refit with `control = sdmTMBcontrol(newton_loops = 1)`
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
m_1 <- run_extra_optimization(m, nlminb_loops = 0, newton_loops = 1)
sanity(m_1)
#> ✓ Non-linear minimizer suggests successful convergence
#> ✓ Hessian matrix is positive definite
#> ✓ No extreme or very small eigen values detected
#> ✓ No gradients with respect to fixed effects are >= 0.001
#> ✓ No fixed-effect standard errors are NA
#> ✓ No fixed-effect standard errors look unreasonably large
#> ✓ No sigma parameters are < 0.01
#> ✓ No sigma parameters are > 100
#> ✓ Range parameter doesn't look unreasonably large
mcmc_res <- residuals(m_1, type = "mle-mcmc", mcmc_iter = 5000, mcmc_warmup = 2000,
stan_args = list(control = list(max_treedepth = 11)))
#> Warning: There were 5 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 11. See
#> http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
#> Warning: Examine the pairs() plot to diagnose sampling problems
#> Warning: The largest R-hat is 1.06, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
png(file = "figures/supp/density/qq.png", units = "in", width = 6.5, height = 6.5, res = 300)
qqnorm(mcmc_res);qqline(mcmc_res)
dev.off()
d$residuals <- mcmc_res[, 1]
# Residuals on map
plot_map_raster +
geom_point(data = d, aes(x = X * 1000, y = Y * 1000, color = residuals)) +
scale_colour_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme(strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90"))
#> Warning: Removed 8 rows containing missing values (geom_point).
ggsave("figures/supp/density/spatial_resid.png", width = 6.5, height = 6.5, dpi = 600)
#> Warning: Removed 8 rows containing missing values (geom_point).
rho is ar_phi on the -1 to 1 scale):tidy(m_1, effects = "ran_pars", conf.int = TRUE)
predict_mcod <- predict(m_1, newdata = pred_grid)
# Save prediction as csv (smaller than the model object...)
write.csv(predict_mcod, "output/predict_mcod_density.csv")
# Plot predicted density and random effects
plot_map_raster +
geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 5) +
labs(fill = expression(kg/km^2)) +
ggtitle("Predicted density (fixed + random)")
#> filter: removed 36,936 rows (15%), 212,517 rows remaining
# I save this below instead to better match the main figure
#ggsave("figures/supp/density/est_map.png", width = 6.5, height = 6.5, dpi = 600)
# Plot spatiotemporal random effect
plot_map_raster +
geom_raster(data = predict_mcod, aes(x = X * 1000, y = Y * 1000, fill = epsilon_st)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
ggtitle("Spatiotemporal random effects")
ggsave("figures/supp/density/epsilon_st_map.png", width = 6.5, height = 6.5, dpi = 600)
# Plot spatial random effect
plot_map_raster +
geom_raster(data = filter(predict_mcod, year == 1999), aes(x = X * 1000, y = Y * 1000, fill = omega_s)) +
scale_fill_gradient2() +
facet_wrap(~ year, ncol = 5) +
theme_plot() +
ggtitle("Spatial random effects")
#> filter: removed 240,214 rows (96%), 9,239 rows remaining
ggsave("figures/supp/density/omega_s_map.png", width = 6.5, height = 6.5, dpi = 600)
get_index_sims method to avoid slow bias correction (not done above!)# preds_cod_sim <- predict(m_1, newdata = pred_grid, sims = 100)
# x <- get_index_sims(preds_cod_ave_sim, area = rep(2*2, nrow(preds_cod_ave_sim)))
#
# ggplot(x, aes(year, est/1000, ymin = lwr/1000, ymax = upr/1000)) +
# geom_line() +
# geom_ribbon(alpha = 0.4)
# x_sims <- get_index_sims(preds_cod_ave_sim, return_sims = TRUE)
# ggplot(x_sims, aes(as.factor(year), log(.value))) +
# geom_violin()
###
# Now get the total index by specifying the area of a grid cell
# Total prediction
preds_cod_sim <- predict(m_1, newdata = pred_grid, sims = 100)
#> Warning: The `sims` argument of `predict.sdmTMB()` is deprecated as of sdmTMB 0.0.21.
#> Please use the `nsim` argument instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
# Total prediction (omitting SD 24!)
preds_cod25_28_sim <- predict(m_1, newdata = filter(pred_grid, !sub_div == 24), sims = 100)
#> filter: removed 38,421 rows (15%), 211,032 rows remaining
# SD 24-25
preds_cod24_25_sim <- predict(m_1, newdata = filter(pred_grid, sub_div %in% c(24, 25)), sims = 100)
#> filter: removed 140,562 rows (56%), 108,891 rows remaining
# SD 24
preds_cod24_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 24), sims = 100)
#> filter: removed 211,032 rows (85%), 38,421 rows remaining
# SD 25
preds_cod25_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 25), sims = 100)
#> filter: removed 178,983 rows (72%), 70,470 rows remaining
# SD 26
preds_cod26_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 26), sims = 100)
#> filter: removed 184,410 rows (74%), 65,043 rows remaining
# SD 27
preds_cod27_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 27), sims = 100)
#> filter: removed 224,505 rows (90%), 24,948 rows remaining
# SD 28
preds_cod28_sim <- predict(m_1, newdata = filter(pred_grid, sub_div == 28), sims = 100)
#> filter: removed 198,882 rows (80%), 50,571 rows remaining
# Now calculate the CPUE index (average)
index_sim <- get_index_sims(preds_cod_sim, area = rep(2*2, nrow(preds_cod_sim)))
index24_25_sim <- get_index_sims(preds_cod24_25_sim, area = rep(2*2, nrow(preds_cod24_25_sim)))
index25_28_sim <- get_index_sims(preds_cod25_28_sim, area = rep(2*2, nrow(preds_cod25_28_sim)))
index24_sim <- get_index_sims(preds_cod24_sim, area = rep(2*2, nrow(preds_cod24_sim)))
index25_sim <- get_index_sims(preds_cod25_sim, area = rep(2*2, nrow(preds_cod25_sim)))
index26_sim <- get_index_sims(preds_cod26_sim, area = rep(2*2, nrow(preds_cod26_sim)))
index27_sim <- get_index_sims(preds_cod27_sim, area = rep(2*2, nrow(preds_cod27_sim)))
index28_sim <- get_index_sims(preds_cod28_sim, area = rep(2*2, nrow(preds_cod28_sim)))
# Join indices
index_sim <- index_sim %>% mutate(sub_div = "Total")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index24_sim <- index24_sim %>% mutate(sub_div = "24")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index25_sim <- index25_sim %>% mutate(sub_div = "25")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index26_sim <- index26_sim %>% mutate(sub_div = "26")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index27_sim <- index27_sim %>% mutate(sub_div = "27")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
index28_sim <- index28_sim %>% mutate(sub_div = "28")
#> mutate: new variable 'sub_div' (character) with one unique value and 0% NA
div_index_sim <- bind_rows(index_sim, index24_sim, index25_sim, index26_sim, index27_sim, index28_sim)
# Now do the same but excluding the areas not sampled in the pred grid
ncells_sub <- filter(pred_grid, year == max(pred_grid$year) & depth > 20 & depth < 120) %>% nrow()
#> filter: removed 242,779 rows (97%), 6,674 rows remaining
pred_grid_sub <- filter(pred_grid, depth > 20 & depth < 120)
#> filter: removed 69,255 rows (28%), 180,198 rows remaining
pred_avg_sim_sub <- predict(m_1, newdata = pred_grid_sub, nsim = 100)
index_avg_sim_sub <- get_index_sims(pred_avg_sim_sub, area = rep(2*2, nrow(pred_avg_sim_sub)))
# Combine and plot & compare
index_avg_sim_comp <- bind_rows(mutate(index_sim, area = "full"),
mutate(index_avg_sim_sub, area = "subset")) %>%
mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> mutate: new variable 'est_t' (double) with 54 unique values and 0% NA
#> new variable 'lwr_t' (double) with 54 unique values and 0% NA
#> new variable 'upr_t' (double) with 54 unique values and 0% NA
ggplot(index_avg_sim_comp, aes(year, est_t, ymin = lwr_t, ymax = upr_t, color = area)) +
geom_line() +
geom_ribbon(alpha = 0.4)
div_index_sim <- div_index_sim %>% mutate(est_t = est/1000, lwr_t = lwr/1000, upr_t = upr/1000)
ind_plot_sd <-
ggplot() +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
geom_line(data = filter(div_index_sim, !sub_div == "Total"), aes(year, est_t, color = sub_div)) +
geom_ribbon(data = filter(div_index_sim, !sub_div == "Total"), aes(year, ymin = lwr_t, ymax = upr_t, fill = sub_div), alpha = 0.2) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2", name = "ICES Sub division") +
guides(color = "none") +
labs(x = "Year", y = "Tonnes", fill = "Sub Divisions") +
theme_plot() +
theme(axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.25, 0.15, 0), "cm")) +
NULL
ind_plot_tot <-
ggplot() +
geom_line(data = filter(div_index_sim, sub_div == "Total"),
aes(year, est_t)) +
geom_ribbon(data = filter(div_index_sim, sub_div == "Total"),
aes(year, ymin = lwr_t, ymax = upr_t, fill = "Total"), alpha = 0.4) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_fill_manual(values = "grey40") +
labs(x = "Year", y = "Tonnes", fill = "") +
theme_plot() +
theme(axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.key.width = unit(0.2, "cm"),
legend.background = element_blank(),
plot.margin = unit(c(0, 0.2, 0.15, 0), "cm")) +
NULL
cpue_map_94_18 <- plot_map_raster_labels +
geom_raster(data = filter(predict_mcod, year %in% c("1994", "2018") & depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 2) +
labs(fill = expression(kg/km^2)) +
theme_plot() +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.text = element_text(size = 5, angle = 90),
legend.position = "bottom") +
NULL
(ind_plot_sd | ind_plot_tot) / cpue_map_94_18 + plot_annotation(tag_levels = 'A')
ggsave("figures/density.png", width = 6.5, height = 6.5, dpi = 600)
# All years cpue
plot_map_raster +
geom_raster(data = filter(predict_mcod, depth < 120 & depth > 10), aes(x = X * 1000, y = Y * 1000, fill = exp(est))) +
scale_fill_viridis_c(trans = "sqrt") +
facet_wrap(~ year, ncol = 5) +
labs(fill = expression(kg/km^2)) +
theme_plot() +
theme(plot.margin = unit(c(0, 0, 0, 0), "cm"),
legend.text = element_text(size = 5, angle = 90),
legend.position = "bottom") +
NULL
ggsave("figures/supp/density/all_years_condition_map_covars.png", width = 6.5, height = 6.5, dpi = 600)
# Temporal trends of sprat by sd
pred_grid %>%
ggplot(aes(x = year, y = biomass_spr_sd, color = factor(sub_div))) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
geom_point(size = 2.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2", name = "Ices sub division") +
guides(fill = "none", linetype = "none") +
labs(y = "Sprat biomass [tonnes]", x = "Year",
color = "") +
theme_plot() +
#facet_wrap(~sub_div) +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = c(0.8, 0.8),
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0.4, 0, 0), "cm")) +
NULL
#> Warning: Removed 4458 rows containing non-finite values (stat_smooth).
#> Warning: Removed 4458 rows containing missing values (geom_point).
ggsave("figures/supp/sprat_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
#> Warning: Removed 4458 rows containing non-finite values (stat_smooth).
#> Removed 4458 rows containing missing values (geom_point).
# Plot bathymetry
bathy_plot <- plot_map_raster_labels +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = depth)) +
scale_fill_viridis(option = "A", direction = -1) +
labs(fill = "Depth") +
theme_plot() +
theme(legend.position = "bottom",
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
#> filter: removed 240,214 rows (96%), 9,239 rows remaining
# Now calculate quantiles of the biomass-weighted values
pal <- c(rev(brewer.pal(n = 8, name = "Dark2")), "gray20")
wm_depth <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.75))) %>%
pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
names_to = "series", values_to = "depth") %>%
group_by(series) %>% # standardize within for easy plotting
mutate(depth_ct = depth - mean(depth)) %>%
ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st quartile, 3rd quartile) into (series, depth) [was 27x4, now 81x3]
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'depth_ct' (double) with 72 unique values and 0% NA
#> ungroup: no grouping variables
plot_w_dep <- ggplot(wm_depth, aes(year, depth, color = series, group = series, fill = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
geom_point(size = 1.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
labs(y = "Depth [m]", x = "Year", color = "") +
scale_y_reverse() +
theme_plot() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0.4, 0, 0), "cm")) +
NULL
# Plot oxygen in 2006
oxy_plot <- plot_map_raster_labels +
geom_raster(data = filter(predict_mcod, year == "2006"), aes(x = X * 1000, y = Y * 1000, fill = oxy)) +
scale_fill_viridis() +
labs(fill = expression(paste("O" [2], " [ml/L]", sep = ""))) +
theme_plot() +
theme(legend.position = "bottom",
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
#> filter: removed 240,214 rows (96%), 9,239 rows remaining
# Plot biomass-weighted oxygen
# Weighted total
wm_oxy <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
names_to = "series", values_to = "oxy") %>%
ungroup() %>%
group_by(series) %>% # standardize within for easy plotting
mutate(oxy_ct = oxy - mean(oxy)) %>%
ungroup()
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 4 columns, ungrouped
#> pivot_longer: reorganized (Median, 1st quartile, 3rd quartile) into (series, oxy) [was 27x4, now 81x3]
#> ungroup: no grouping variables
#> group_by: one grouping variable (series)
#> mutate (grouped): new variable 'oxy_ct' (double) with 81 unique values and 0% NA
#> ungroup: no grouping variables
# Find which depths to calculate average oxygen on
wm_depth %>% filter(series == "1st quartile") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
wm_depth %>% filter(series == "3rd quartile") %>% summarise(mean_depth = mean(depth))
#> filter: removed 54 rows (67%), 27 rows remaining
#> summarise: now one row and one column, ungrouped
annual_oxy <- pred_grid %>%
drop_na(oxy) %>%
filter(depth > 29 & depth < 61) %>% ###
#filter(depth > 40 & depth < 60) %>% ###
group_by(year) %>%
summarize(median_oxy = median(oxy)) %>%
ungroup() %>%
rename("oxy" = "median_oxy") %>%
mutate(series = "Median (env.)",
oxy_ct = oxy - mean(oxy))
#> drop_na: no rows removed
#> filter: removed 179,847 rows (72%), 69,606 rows remaining
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> rename: renamed one variable (oxy)
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'oxy_ct' (double) with 27 unique values and 0% NA
oxygen_series <- bind_rows(wm_oxy, annual_oxy)
plot_w_oxy <-
ggplot(oxygen_series, aes(year, oxy, color = series, group = series, fill = series, linetype = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
geom_point(size = 1.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_linetype_manual(values = c(1,1,1,2)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
color = "") +
theme_plot() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
plot_w_oxy
(bathy_plot | oxy_plot) / (plot_w_dep | plot_w_oxy) +
plot_annotation(tag_levels = 'A')
ggsave("figures/weighted_quantiles.png", width = 6.5, height = 6.5, dpi = 600)
# Now see how the weighted mean oxygen differs from the average oxygen at the range of depths used by cod
wm_depth %>% group_by(series) %>% summarise(mean_depth = mean(depth))
#> group_by: one grouping variable (series)
#> summarise: now 3 rows and 2 columns, ungrouped
wm_oxy2 <- wm_oxy %>% filter(series == "Median") %>% mutate(series = "Biomass-weighted")
#> filter: removed 54 rows (67%), 27 rows remaining
#> mutate: changed 27 values (100%) of 'series' (0 new NA)
pred_grid_ave_oxy <- pred_grid %>%
filter(depth < 60 & depth > 30) %>%
group_by(year) %>%
summarise(oxy = mean(oxy)) %>%
mutate(series = "Mean oxygen at depth 30-60 m")
#> filter: removed 185,247 rows (74%), 64,206 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> mutate: new variable 'series' (character) with one unique value and 0% NA
combi <- bind_rows(wm_oxy2, pred_grid_ave_oxy)
ggplot(combi, aes(year, oxy, color = series)) +
geom_point() +
labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
color = "") +
scale_color_manual(values = pal[c(1,3)], name = "") +
stat_smooth(method = "lm") +
theme_plot()
#> `geom_smooth()` using formula 'y ~ x'
ggsave("figures/supp/density/weighted_oxy_vs_depth_oxy.png", width = 6.5, height = 6.5, dpi = 600)
#> `geom_smooth()` using formula 'y ~ x'
# Oxygen
wm_oxy_sd <- predict_mcod %>%
group_by(year, sub_div) %>%
summarise("weighted_median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)))
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 3 columns, one group variable remaining (year)
ggplot(wm_oxy_sd, aes(year, weighted_median, color = factor(sub_div))) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
geom_point(size = 2.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
guides(color = "none", fill = "none", linetype = "none") +
labs(y = expression(paste("O" [2], " [ml/L]", sep = "")), x = "Year",
color = "") +
theme_plot() +
facet_wrap(~sub_div) +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
ggsave("figures/supp/density/weighted_oxygen_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
# # Calculate weighted mean and quantiles for sd 25 only
predict_mcod %>%
filter(sub_div == 25 & year >= 2015) %>%
#group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.25)),
"3rd quartile" = hutils::weighted_quantile(v = oxy, w = exp(est), p = c(0.75))) %>%
as.data.frame()
#> filter: removed 236,403 rows (95%), 13,050 rows remaining
#> summarise: now one row and 3 columns, ungrouped
# Now compare delta change in oxygen with delta change in condition
cond_delta_sd <- read.csv("output/div_index_delta_cond.csv") %>%
mutate(sub_div = as.character(sub_div)) %>%
dplyr::select(-X) %>%
filter(!sub_div == "Total")
#> mutate: no changes
#> filter: removed one row (17%), 5 rows remaining
wm_oxy_sd_delta <- wm_oxy_sd %>%
group_by(sub_div) %>%
filter(year %in% c(1993, 1994, 1995, 2017, 2018, 2019)) %>%
dplyr::select(year, weighted_median, sub_div) %>%
pivot_wider(1:3, names_from = year, values_from = weighted_median) %>%
mutate(start_oxy = mean(`1993`, `1994`, `1995`),
end_oxy = mean(`2017`, `2018`, `2019`),
delta_oxy = end_oxy - start_oxy) %>%
ungroup() %>%
dplyr::select(sub_div, delta_oxy) %>%
mutate(sub_div = as.character(sub_div))
#> group_by: one grouping variable (sub_div)
#> filter (grouped): removed 105 rows (78%), 30 rows remaining
#> pivot_wider: reorganized (year, weighted_median) into (1993, 1994, 1995, 2017, 2018, …) [was 30x3, now 5x7]
#> mutate (grouped): new variable 'start_oxy' (double) with 5 unique values and 0% NA
#> new variable 'end_oxy' (double) with 5 unique values and 0% NA
#> new variable 'delta_oxy' (double) with 5 unique values and 0% NA
#> ungroup: no grouping variables
#> mutate: converted 'sub_div' from double to character (0 new NA)
# delta_dat <- cond_delta_sd %>% left_join(wm_oxy_sd_delta)
#
# ggplot(delta_dat, aes(delta_oxy, delta_cond, color = sub_div)) +
# geom_point(size = 3) +
# geom_vline(xintercept = 0, color = "red", linetype = 2) +
# geom_hline(yintercept = 0, color = "red", linetype = 2)
# NULL
# Depth
wm_dep_sd <- predict_mcod %>%
group_by(year, sub_div) %>%
summarise("weighted_median" = hutils::weighted_quantile(v = depth, w = exp(est), p = c(0.5)))
#> group_by: 2 grouping variables (year, sub_div)
#> summarise: now 135 rows and 3 columns, one group variable remaining (year)
ggplot(wm_dep_sd, aes(year, weighted_median, color = factor(sub_div))) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), size = 1) +
geom_point(size = 2.5, alpha = 0.8) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
guides(color = "none", fill = "none", linetype = "none") +
labs(y = "Depth [m]", x = "Year", color = "") +
theme_plot() +
facet_wrap(~sub_div) +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 8),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
ggsave("figures/supp/density/weighted_depth_subdiv.png", width = 6.5, height = 6.5, dpi = 600)
# Plot temperature in 2006
temp_plot <- plot_map_raster +
geom_raster(data = filter(predict_mcod, year == "1999"), aes(x = X * 1000, y = Y * 1000, fill = temp)) +
scale_fill_viridis(option = "inferno", direction = -1) +
labs(fill = "Temperature [°C]") +
theme_plot() +
theme(legend.position = "bottom",
plot.margin = unit(c(0, 0.1, 0, 0.2), "cm")) +
NULL
annual_temp <- pred_grid %>%
drop_na(temp) %>%
group_by(year) %>%
summarize(median_temp = median(temp)) %>%
ungroup() %>%
rename("temp" = "median_temp") %>%
mutate(series = "Median (env.)",
temp_ct = temp - mean(temp))
# Weighted total
wm_temp <- predict_mcod %>%
group_by(year) %>%
summarise("Median" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.5)),
"1st quartile" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.1)),
"3rd quartile" = hutils::weighted_quantile(v = temp, w = exp(est), p = c(0.9))) %>%
pivot_longer(cols = c("Median", "1st quartile", "3rd quartile"),
names_to = "series", values_to = "temp") %>%
ungroup() %>%
group_by(series) %>% # standardize within for easy plotting
mutate(temp_ct = temp - mean(temp)) %>%
ungroup()
temp_series <- bind_rows(wm_temp, annual_temp)
plot_w_temp <- ggplot(temp_series, aes(year, temp, color = series, group = series, fill = series, linetype = series)) +
stat_smooth(method = "gam", formula = y ~ s(x, k = 4), se = FALSE, size = 1) +
geom_point(size = 2.5, alpha = 0.8, color = "white", shape = 21) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 10)) +
scale_linetype_manual(values = c(1,1,1,2)) +
scale_color_manual(values = pal, name = "Quantiles") +
scale_fill_manual(values = pal) +
guides(fill = "none", linetype = "none", color = guide_legend(nrow = 2, title.position = "top", title.hjust = 0.5)) +
labs(y = "Temperature [°C]", x = "Year",
color = "") +
theme_plot() +
theme(aspect.ratio = 1,
axis.text.x = element_text(angle = 0),
legend.position = "bottom",
legend.key.size = unit(0.6, "cm"),
legend.text = element_text(size = 7),
plot.margin = unit(c(0, 0, 0, 0), "cm")) +
NULL
temp_plot + plot_w_temp +
plot_annotation(tag_levels = 'A') + plot_layout(ncol = 2)
ggsave("figures/supp/density/weighted_temp.png", width = 6.5, height = 8.5, dpi = 600)
annual_temp2 <- pred_grid %>%
drop_na(temp) %>%
group_by(year) %>%
summarize(value = mean(temp)) %>%
mutate(series = "Mean (env.)",
variable = "Temperature",
species = "Environment")
#> drop_na: no rows removed
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
annual_oxy2 <- pred_grid %>%
drop_na(oxy) %>%
group_by(year) %>%
summarize(value = mean(oxy)) %>%
ungroup() %>%
mutate(series = "Mean (env.)",
variable = "Oxygen",
species = "Environment")
#> drop_na: no rows removed
#> group_by: one grouping variable (year)
#> summarize: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
# Weighted total
wm_temp2_cod <- predict_mcod %>%
filter(exp(est) < 3000) %>%
group_by(year) %>%
summarise(value = weighted.mean(x = temp, w = exp(est))) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Temperature",
species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
wm_temp2_fle <- predict_mcod %>%
filter(density_fle < 2000) %>%
group_by(year) %>%
summarise(value = weighted.mean(x = temp, w = density_fle)) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Temperature",
species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
wm_temp2_spr <- predict_mcod %>%
group_by(year) %>%
drop_na(biomass_spr) %>%
summarise(value = weighted.mean(x = temp, w = biomass_spr)) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Temperature",
species = "Sprat")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
wm_temp2_her <- predict_mcod %>%
group_by(year) %>%
drop_na(biomass_her) %>%
summarise(value = weighted.mean(x = temp, w = biomass_her)) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Temperature",
species = "Herring")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
wm_oxy2_cod <- predict_mcod %>%
filter(exp(est) < 3000) %>%
group_by(year) %>%
summarise(value = weighted.mean(x = oxy, w = exp(est))) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Oxygen",
species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
wm_oxy2_fle <- predict_mcod %>%
filter(density_fle < 2000) %>%
group_by(year) %>%
summarise(value = weighted.mean(x = oxy, w = density_fle)) %>%
ungroup() %>%
mutate(series = "Mean (weighted)",
variable = "Oxygen",
species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 2 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'series' (character) with one unique value and 0% NA
#> new variable 'variable' (character) with one unique value and 0% NA
#> new variable 'species' (character) with one unique value and 0% NA
mean_oxy_temp <- bind_rows(annual_temp2, annual_oxy2, wm_temp2_cod, wm_temp2_fle, wm_oxy2_cod, wm_oxy2_fle, wm_temp2_her, wm_temp2_spr)
# coords
wm_coord_cod <- predict_mcod %>%
filter(exp(est) < 3000) %>%
group_by(year) %>%
summarise(wlat = weighted.mean(x = lat, w = exp(est)),
wlon = weighted.mean(x = lon, w = exp(est))) %>%
ungroup() %>%
mutate(species = "Cod")
#> filter: removed 714 rows (<1%), 248,739 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA
wm_coord_fle <- predict_mcod %>%
filter(density_fle < 2000) %>%
group_by(year) %>%
summarise(wlat = weighted.mean(x = lat, w = density_fle),
wlon = weighted.mean(x = lon, w = density_fle)) %>%
ungroup() %>%
mutate(species = "Flounder")
#> filter: removed 1,291 rows (1%), 248,162 rows remaining
#> group_by: one grouping variable (year)
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA
wm_coord_spr <- predict_mcod %>%
group_by(year) %>%
drop_na(biomass_spr) %>%
summarise(wlat = weighted.mean(x = lat, w = biomass_spr),
wlon = weighted.mean(x = lon, w = biomass_spr)) %>%
ungroup() %>%
mutate(species = "Sprat")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA
wm_coord_her <- predict_mcod %>%
group_by(year) %>%
drop_na(biomass_her) %>%
summarise(wlat = weighted.mean(x = lat, w = biomass_her),
wlon = weighted.mean(x = lon, w = biomass_her)) %>%
ungroup() %>%
mutate(species = "Herring")
#> group_by: one grouping variable (year)
#> drop_na (grouped): removed 21,502 rows (9%), 227,951 rows remaining
#> summarise: now 27 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'species' (character) with one unique value and 0% NA
coord_data <- bind_rows(wm_coord_cod, wm_coord_fle, wm_coord_spr, wm_coord_her)
# Plot
palz <- brewer.pal(n = 5, name = "Dark2")
palz[2] <- "gray50"
mean_oxy_temp %>%
filter(variable == "Temperature",
!year == 1994) %>%
ggplot(aes(year, value, color = species, linetype = species, size = species)) +
geom_point(size = 2.3, alpha = 0.5) +
guides(size = "none", linetype = "none") +
scale_linetype_manual(values = c(1,2,1,1,1)) +
scale_size_manual(values = c(1,2.3,1,1,1)) +
scale_color_manual(values = palz, name = "") +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3), se = FALSE) +
labs(x = "Year", y = "Temperature [°C]") +
theme(legend.position = "bottom",
aspect.ratio = 1) +
NULL
#> filter: removed 86 rows (40%), 130 rows remaining
ggsave("figures/supp/density/weighted_mean_temp_species.png", width = 6.5, height = 8.5, dpi = 600)
# Now do flounder
ggplot(filter(coord_data, !year == 1993), aes(wlon, wlat, color = year)) +
facet_wrap(~species, scales = "free") +
geom_point() +
geom_path() +
scale_color_viridis()
#> filter: removed 4 rows (4%), 104 rows remaining
coord_data %>% filter(year %in% c(1994, 2019)) %>%
ggplot(aes(wlon, wlat, color = factor(year))) +
facet_wrap(~species, scales = "free") +
geom_point() +
geom_path() +
scale_color_viridis(discrete = TRUE)
#> filter: removed 100 rows (93%), 8 rows remaining
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
#> geom_path: Each group consists of only one observation. Do you need to adjust
#> the group aesthetic?
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
utm_coordz <- LongLatToUTM(zone = 33, coord_data$wlon, coord_data$wlat)
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
#> prefer_proj): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition
coord_data$wlat_y <- utm_coordz$Y
coord_data$wlon_x <- utm_coordz$X
plot_map_raster +
facet_wrap(~species) +
geom_point(data = filter(coord_data, year %in% c(1994, 2019)),
aes(wlon_x, wlat_y, color = factor(year)), inherit.aes = FALSE, size = 3) +
theme_plot() +
scale_color_brewer(palette = "Dark2", name = "Year") +
theme(legend.position = "bottom",
plot.margin = unit(c(0, 0.1, 0, 0.2), "cm")) +
NULL
#> filter: removed 100 rows (93%), 8 rows remaining
plot_map_raster +
facet_wrap(~species) +
geom_point(data = filter(coord_data, !year == c(1993)),
aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 1) +
geom_path(data = filter(coord_data, !year == c(1993)),
aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 0.1) +
theme_plot() +
scale_color_viridis() +
theme(legend.position = "bottom",
plot.margin = unit(c(0, 0.1, 0, 0.2), "cm"),
legend.text = element_text(angle = 30)) +
xlim(1.8*xmin2, 0.85*xmax2) +
ylim(1.023*ymin2, 0.976*ymax2) +
NULL
#> filter: removed 4 rows (4%), 104 rows remaining
#> filter: removed 4 rows (4%), 104 rows remaining
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
#> Scale for 'y' is already present. Adding another scale for 'y', which will
#> replace the existing scale.
ggsave("figures/supp/density/weighted_mean_coord_species.png", width = 6.5, height = 8.5, dpi = 600)
# pred_grid %>%
# filter(year > 2010 & sub_div == 25 & depth > 48 & depth < 52) %>%
# group_by(year) %>%
# summarise(mean_oxy = mean(oxy))
#
#
# pred_grid %>%
# filter(year > 2010 & sub_div == 25 & depth > 50 & depth < 65) %>%
# group_by(year) %>%
# summarise(mean_oxy = mean(oxy))
# VR combined plot
p1 <- plot_map_raster +
facet_wrap(~species) +
geom_point(data = filter(coord_data, !year == c(1993)),
aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 1) +
geom_path(data = filter(coord_data, !year == c(1993)),
aes(wlon_x, wlat_y, color = year), inherit.aes = FALSE, size = 0.1) +
theme_plot() +
scale_color_viridis() +
theme(legend.position = "bottom",
legend.spacing.x = unit(0, 'cm'),
legend.spacing.y = unit(0, 'cm'),
legend.margin=margin(0,0, 0, 0),
legend.box.margin=margin(-20,-20,-20,-20),
#strip.placement = NULL,
#plot.margin = unit(c(0, 0.1, 0, 0.2), "cm"),
legend.text = element_text(angle = 30)#,
#plot.margin = unit(c(1, 1, 1, 1), "cm")
) +
xlim(1.8*xmin2, 0.85*xmax2) +
ylim(1.023*ymin2, 0.976*ymax2) +
NULL
#> filter: removed 4 rows (4%), 104 rows remaining
#> filter: removed 4 rows (4%), 104 rows remaining
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
#> Scale for 'y' is already present. Adding another scale for 'y', which will
#> replace the existing scale.
p2 <- mean_oxy_temp %>%
filter(variable == "Temperature",
!year == 1994) %>%
ggplot(aes(year, value, color = species, linetype = species, size = species)) +
geom_point(size = 2.3, alpha = 0.5) +
guides(size = "none", linetype = "none", color = guide_legend(nrow = 1)) +
scale_linetype_manual(values = c(1,2,1,1,1)) +
scale_size_manual(values = c(1,2.3,1,1,1)) +
scale_color_manual(values = palz, name = "") +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3), se = FALSE) +
labs(x = "Year", y = "Temperature [°C]") +
theme(legend.position = "bottom",
legend.text = element_text(size = 6),
legend.spacing.x = unit(0, 'cm'),
legend.spacing.y = unit(0, 'cm'),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-40,-40,-40,-40),
aspect.ratio = 1#,
#plot.margin = unit(c(1, 1, 1, 1), "cm")
) +
NULL
#> filter: removed 86 rows (40%), 130 rows remaining
p1 | p2
ggsave("figures/supp/density/VR_weighted_mean_coord_species.png", width = 6.5, height = 8.5, dpi = 600)
# Prepare prediction data frame
d2 <- d %>% drop_na(oxy)
nd_oxy <- data.frame(oxy = seq(min(d2$oxy), max(d2$oxy), length.out = 100))
nd_oxy <- nd_oxy %>%
mutate(year = 2003L,
depth_sc = 0,
oxy_sc = (oxy - mean(oxy))/sd(oxy),
temp_sc = 0)
# Predict
p_margin_oxy <- predict(m_1, newdata = nd_oxy, se_fit = TRUE, re_form = NA)
mar_oxy <- ggplot(p_margin_oxy, aes(oxy, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Oxygen[haul]))
# Prepare prediction data frame
nd_dep <- data.frame(depth = seq(min(d2$depth), max(d2$depth), length.out = 100))
nd_dep <- nd_dep %>%
mutate(year = 2003L,
depth_sc = (depth - mean(depth))/sd(depth),
oxy_sc = 0,
temp_sc = 0)
# Predict
p_margin_dep <- predict(m_1, newdata = nd_dep, se_fit = TRUE, re_form = NA)
mar_depth <- ggplot(p_margin_dep, aes(depth, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Depth[haul]))
# Prepare prediction data frame
nd_temp <- data.frame(temp = seq(min(d$temp), max(d$temp), length.out = 100))
nd_temp <- nd_temp %>%
mutate(year = 2003L,
depth_sc = 0,
oxy_sc = 0,
temp_sc = (temp - mean(temp))/sd(temp))
# Predict
p_margin_temp <- predict(m_1, newdata = nd_temp, se_fit = TRUE, re_form = NA)
mar_temp <- ggplot(p_margin_temp, aes(temp, exp(est),
ymin = exp(est) - 1.96 * exp(est_se), ymax = exp(est) + 1.96 * exp(est_se))) +
geom_line() +
geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(30, 220)) +
theme(aspect.ratio = 1) +
labs(x = expression(Temp[haul]))
# mar_temp + mar_depth + mar_oxy +
# plot_annotation(tag_levels = 'A')
#
# ggsave("figures/supp/density/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)
p_margin_dep2 <- p_margin_dep %>%
mutate(variable = "Depth") %>%
dplyr::select(est, est_se, depth_sc, variable) %>%
rename(value = depth_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_oxy2 <- p_margin_oxy %>%
mutate(variable = "Oxygen") %>%
dplyr::select(est, est_se, oxy_sc, variable) %>%
rename(value = oxy_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
p_margin_tem2 <- p_margin_temp %>%
mutate(variable = "Temperature") %>%
dplyr::select(est, est_se, temp_sc, variable) %>%
rename(value = temp_sc)
#> mutate: new variable 'variable' (character) with one unique value and 0% NA
#> rename: renamed one variable (value)
margs <- bind_rows(p_margin_dep2, p_margin_oxy2, p_margin_tem2)
ggplot(margs, aes(value, exp(est), ymin = exp(est - 1.96 * est_se),
ymax = exp(est + 1.96 * est_se))) +
geom_ribbon(alpha = 0.4) +
geom_line() +
facet_wrap(~variable, scales = "free_x") +
labs(y = "Biomass density",
x = "Scaled value") +
theme_plot() +
theme(axis.text.x = element_text(angle = 0)) +
NULL
ggsave("figures/supp/density/marginal_effects.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()